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An algorithm is presented which computes a translationally invariant matrix product state ap- 
proximation of the ground state of an infinite ID system; it does this by embedding sites into 
an approximation of the infinite "environment" of the chain, allowing the sites to relax, and then 
merging them with the environment in order to refine the approximation. By making use of matrix 
product operators, our approach is able to directly model any long-range interaction that can be 
systematically approximated by a series of decaying exponentials. We apply these techniques to 
compute the ground state of the Haldane-Shastry model and present results. 
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I. INTRODUCTION 

Typical macroscopic quantities of material contain par- 
ticules numbering on the order of 10 26 . The only hope 
of modeling such systems is to come up with a repre- 
sentation that is relatively small but nonetheless able to 
capture essential properties. For quantum systems this 
is particularly difficult because the state space grows ex- 
ponentially with the size of the system, so there is a large 
amount of information which needs to be disposed. 

One approach is to start with a small system and then 
build it up, using a renormalization transformation at 
each step to project out "undesirable" states and keep 
the representation small. The key is to get one's crite- 
ria for "undesirable" correct. When studying low-energy 
behavior of a system, it turns out that the best selec- 
tion criteria is not the energy (unless one is solving the 
Kondo problem^), but the weight in a density matrix, 



which incorporates interactions with a surrounding envi- 
ronment. This idea, first proposed by White^ in his den- 
sity matrix renormalization group (DMRG) algorithm, 
has proven very successful in building effective but com- 
pact representations of large one-dimensional (ID) sys- 
tems. (See Ref. 0for an excellent review of this subject.) 
However, it has not worked as well as desired in mod- 
eling systems with long-range interactions, and it does 
not produce manifestly translationally invariant repre- 
sentations of translationally invariant states. One can 
improve on both these fronts by working in momentum 
space rather than in real spaced, but the computation 
is more costly as the transformation to momentum space 
introduces many additional operators that need to be 
renormalized at each step, and it can break symmetries 
that had been present in real space. Also, none of these 
methods are particularly effective in two or more dimen- 
sions. 

An alternative approach is inspired by observing that 
the DMRG technique converges to a matrix product state 
(MPS). 6 Rather than building up progressively larger 
systems, the iTEBD algorithm 7 assumes a priori that the 
system takes an infinite translationally invariant MPS 
form, and then uses imaginary time evolution to con- 
verge to the ground state. This approach has the ad- 
vantage that it obtains a very compact representation - 
in particular, it dispenses with the need to form renor- 
malized versions of operators. Furthermore, it admits a 
natural extension to two dimensions 8 by using projected 
entangled-pair states (PEPS)^, a two-dimensional gen- 
eralization of MPS. However, it is only able to handle 
systems with short-range interactions. 

In this paper, we present an algorithm which also com- 
putes a translationally invariant MPS representation of 
a ground state, but without being limited to short-range 
interactions. Instead of simulating an evolution in imag- 
inary time as in the iTEBD algorithm, our approach 
follows the spirit of the variational technique described 
in Ref. [l(| but enhanced to apply to infinite systems; 
additionally, we build in the ability to work with arbi- 
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trary matrix product operator s 11 1 12 , which makes this al- 
gorithm naturally suited to handle systems with long- 
range interactions. 

The intuition behind our approach is as follows. Sup- 
pose that one were given an infinitely large, translation- 
ally invariant ID quantum chain held at zero tempera- 
ture. If one were to add an additional site to this chain 
and allow the chain to relax, then one would expect that 
all of the old sites would remain unchanged, and the new 
site would change to match the rest. If one could emu- 
late the environment experienced by a single site in this 
infinite chain, then by embedding a site into this environ- 
ment and allowing the system to relax, one would obtain 
a site that "looks like" all of the sites in our infinite chain, 
giving us a compact representation for the chain. 

We note that this algorithm bears some similarity 
to the product wave function renormalization group 
(PWFRG) 13 in that both algorithms have the same 
goal — to compute representations of ground states for in- 
finite systems — and both use the same underlying matrix 
product structure to represent states^. The difference is 
that while the PWFRG approach seeks the infinite limit 
by starting with a small system and progressively enlarg- 
ing it, we start from the very beginning with an infinite 
system, represented in terms of effective environments 
which we progressively refine. Furthermore, we incorpo- 
rate matrix product operators into our approach which 
allow us to model systems with long-range interactions. 

The remainder of this paper shall be organized as fol- 
lows. First, we shall review the matrix product formal- 
ism, and show how it allows us to construct a representa- 
tion of a site embedded in and entangled with an environ- 
ment. Second, we shall outline the operation of an algo- 
rithm which uses an iteration procedure to approximate 
the effective environment of a site in an infinite chain; 
we shall then explain how one can obtain the expected 
values of operators from the output of this algorithm. 
Third, we shall explain how to obtain matrix product 
factorizations of Hamiltonians with exponentially decay- 
ing interactions. Finally, we shall pull all of these ideas 
together and show how they can be applied to obtain an 
accurate MPS representation of the ground state of the 
Haldane-Shastry model — an exactly solvable model with 
long range interactions. 

II. ALGORITHMS 

A. Finding a MPS representation of the ground 
state 

We start by recalling the form of a matrix product 
state for a finite system^. For a system with n sites, a 
matrix product state with boundaries takes the form 

Zl ( L o)i 1 ( S l)?X ( S z)i 2 % ■ ■ ■ ( R n+l) lrl + 1 



The left-hand-side corresponds to the rank-n tensor de- 
scribing our state; each index ak corresponds to a basis of 
a physical observable, such as the z-component of spin, at 
site k. The rank-3 tensors Sk are the site tensors, which 
have two kinds of indices: the superscript index, which 
has dimension d and corresponds to the actual physical 
observable, and the subscript indices, which have dimen- 
sion x an d give information about entanglement between 
each site and its neighbors. The vectors Lq and R n +i give 
the boundary conditions. 

The advantage of the matrix product representation is 
that it decomposes the quantum state into a collection 
of tensors, each of which is naturally associated with a 
site on the lattice. Because of this, we shall see that 
we can "zoom in" on one site and then we only have to 
work with its corresponding tensor and a relatively small 
environment. In this respect, it is almost as nice as a 
truly local description that would let us ignore all the 
other sites entirely-a fact which is remarkable given that 
it can still capture non-local properties resulting from 
entanglement. 

Given this form, we now consider how to compute the 
expected values of operators. As discussed in Ref. 11 and 
Ref. [HI, operators of interest can typically be expressed 
as a matrix product operator, 

Q(a 1 ,a' 1 ),(a 2 ,a' 2 ),...,( a n,a' n ) _ 
{ik} 

Having decomposed both the state and the operator into 
a product of tensors associated with lattice sites, we like- 
wise decompose the expected value into a product of ten- 
sors by combining the state and operator tensor at each 
site to define, 

(eP) =(eP) 

V K Ji,j \ J (i',i",i">), {j'.j" ,r") 

a k ,a' k 

(Note that we have taken all of the left and right subscript 
indices and grouped them together.) To complete our 
decomposition of the expected value, we shall also need 
to define the left and right boundaries, 

K 1 ),.,,,,,,,, =- m, it)* ■ » 

('^'+ 1 )j={ji j» 3 w) := C^n+Oj" (-^n+l)j»/ • 

These boundary vectors can be thought of as defining an 
environment into which our system of n sites has been 
embedded. With these tensors so defined, the expected 
value of the matrix product operator O with respect to 
the matrix product state S is given by (S \ 0\ S) = LP • 
E[ 0) ■ EP ■ ■ ■ E P ■ RP } , where • indicates summation 
over the adjoining subscript indices. 
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Now we have almost arrived at the picture described 
in the introduction, except that we have an environment 
and multiple sites rather than a single site. To compute 
the effective environment of a particular site in this chain, 

we absorb all of the sites surrounding it into the system 

(o) 

environment by contracting the E k matrices into the 
left and right boundaries - that is, we make the inductive 



definitions hf" 1 :- 



Lf_{.EP and R k := E p.Ri% 



By doing this, we can write the expected value of the 
operator as a function of this site tensor, 

{S\G\S)(S k ) = S* k oM^*oS k := (3) 



E 



where o denotes summation over the appropriate adja- 
cent subscript and superscript indices, and 

(«),(«') 



E(4 o, ,) (i ,, 1 „, 1 ,„ ) (0 t )^(4?, 



U',j",j'") 



We have now obtained an explicit means to compute 
the expected value of an observable O as a function of the 
site tensor at position k knowing only the environment of 

(o) (o) 

the site as given by L k _ x and R k+1 - However, recall that 
we want to be more than passive observers — we want to 
actively move the system as close as possible to its ground 
state. Thus, we now want to vary the site tensor Sk in 
order to minimize the energy of the system. If we let the 
Hamiltonian (matrix product) operator be denoted by Ti. 
and the identity operator by X, then employing equation 
© we see that what we seek is the site tensor which 
minimizes the function 



S** o M^ n h o S k 
Sf a M< z >* o S k ' 



which gives us the (normalized) energy of the system as a 
function of only the site tensor at position k (i.e., assum- 
ing that the environment has been frozen in place) . Since 
this is a Rayleigh quotient, computing the minimizer is 
equivalent to solving a generalized eigenvalue problem. 

There is a subtlety in this procedure, however, which 
is that one needs to take steps to make sure that the 
normalization matrix AT^ 1 ^ is well-conditioned. This 
can be done by imposing the "right-normalization" con- 
dition J2 a ,i ( s i)tj ( s i)ij> = %' on a11 tne sites to 

the left of k, and the "left-normalization" condition 
Ea.j ( s i)ij (Sl)i>j = da' on a11 tne sites to the right of k. 
This ensures that the subscript indices connected to Sk 
are orthonormal, which makes M 1 * 1 ^ well-conditioned. 
We shall discuss how this is done in our algorithm shortly; 
for more information on how the normalization condition 
is used for finite- length systems, see Ref. 

Now we have all of the ingredients that we need to 
build our algorithm: a way of expressing a chain as a site 



1. Set x = 1, and L s = R s = 1. 

2. Compute L\ and R\ for the Hamiltonian (matrix prod- 
uct) operator TL and the identity operator T. (The latter, 
of course, has trivial boundaries.) 

3. Until convergence has been reached: 

(a) Use an eigenvalue solver (such as ARPACK— ) to 
obtain the minimal eigenvalue and corresponding 
eigenvector of the generalized eigenvalue problem 
M {H)k S k = \M {T)k S k . To accelerate convergence, 
feed in Sk-i as a starting estimate for the eigenvec- 
tor. 

(b) If this is an odd- numbered step, then right-normalize 
Sk and contract into the left boundary. Specifically: 

i. Merge the superscript and first subscript index of 
Sk to form a matrix, and compute the singular 
value decomposition (SVD), U ■ E • V. 

ii. Set Sk '■= U ■ V , and ungroup indices to return Sk 
to its original rank-3 shape. 

iii. Compute E k for 7i and T using the normalized 
Sk- 

iv. Contract the site into the left boundary by setting 



L 



Lf>.£f> and<> 



fc+i — - fe - — i-k+i — R i 0) ■ ( This ste P 
"absorbs" the site into the environment.) 



(c) If this is an even-numbered step, perform an analo- 
gous process, but /e/f-normalize Sk by merging the 
superscript and second index together before the 
SVD, and then contract the normalized site into the 
ri(//i£-boundary. 

4. If a better approximation to the ground state is desired, 
then increase x f° r the system and repeat step 3. Increasing 
the dimension of a subscript can be done without altering 
the state by multiplying the adjoining tensors by a x x 
(x + Ax) matrix and its inverse; this allows one to build 
on the work of previous iterations (and hence accelerate 
convergence), rather than having to start from scratch with 
the new x- 



TABLE I: Algorithm to compute a (normalized) transla- 
tionally invariant matrix product state representation of the 
ground state of an infinite chain. 



tensor embedded in an environment, and a way to relax a 
system in this representation constrained so that only the 
site tensor (and not its environment) is changed. How- 
ever, up to this point we have been working with a finite 
system; in order to apply these ideas to infinite systems, 

we shall modify our notation slightly to replace position 

(o) 

labels with iteration labels. That is, we shall let L k and 

(O) 

R k denote the infinite environment at iteration step fc, 
and Sk denote the inserted site at this iteration. With 
this notation, the algorithm to find the ground state of a 
system is given in Table HI It is dominated by the costs 
of absorbing the site and operator tensors into the envi- 
ronment (in step 3b. iv, and in the Lanczos iteration in 
step 3a), which are respectively 0(cdx 3 ) and 0(c 2 d 2 x 2 ), 
with c referring to the auxiliary dimension of the operator 
tensor. 
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We note here that this algorithm is unstable when ap- 
plied directly to systems with anti-ferromagnetic inter- 
actions because ground states of such systems are not 
invariant under translations of one site. Happily, since 
such systems are invariant under translations of two sites, 
there is a simple fix: work with blocks of two spins rather 
than one by setting d = 4 and multiplying two of the 
operator tensors together to form a two-site operator 
tensor^ this only affects the inputs to the algorithm 
and does not require changing the algorithm itself. 



B. Computing expected values of operators 



given by 



The output of the algorithm of Sec. Ill Al is a nor- 
malized site tensor S which gives a translationally in- 
variant representation for the ground state of the sys- 
tem. In order to make this useful, we need to have a 
way to obtain the expected value of operators from it. 
Of course, for extensive observables the expected value 
will be infinite since we have an infinitely large system, 
so instead we seek the more useful quantity of the ex- 
pected value per site. We shall consider two cases of 
operators: local operators and general matrix product 
operators. The basic trick in both cases is to note that 
limjv^oo (E^y ■ v = A N ■ v, where A is a matrix in 
the maximal eigenspace of £<°> - that is, in the infinite 
limit the action of the operator will converge to its action 
on the maximal eigenspace, since its action on all other 
eigenspaces will be negligible by comparison. (Assuming, 
of course, that A - v ^ 0.) 

First we consider local operators, which can be ex- 
pressed as sum of terms of the form eg) 0\ ® O2 <8> 
• • • <S> On <8> J® 00 , such as a magnetic field operator, 
I^igiZ®!® 00 , or a two-point correlator, I®°°®Z(i)I® r ® 
Z <S> 7®°°. Since our system is translationally invariant, 
to compute the expected value per site of this operator 
we need only evaluate the expected value of one term in 
the sum. Assuming that E^ has a non-degenerate max- 
imal eigenvalue, we have that for (almost) any vector v, 
(E^) 00 -v ~ v R and v- (E^) 00 ~ v L , where v L and v R 
are the respective left and right eigenvectors correspond- 
ing to the maximal eigenvalue. Ergo, the expected value 
of one term of this operator (and thus the expected value 
per site) is given by 

v L ■ E Ql ■ E° 2 . . . E° N ■ v R 



(E<?)) 



N 



Next we consider general matrix product operators. 
We start by writing an expression for the expected value 
per site for a finite chain of length N, whose site tensors 
are copies of the (normalized) site tensor obtained from 
the variational algorithm in Sec. Ill Al Since the infinite 
chain has no boundaries, to compute the expected value 
of the finite chain we need to explicitly supply left and 
right boundaries, respectively L s and R s , and so the ex- 
pected value per site of O is a function of the boundaries 



£ N (L S ,R S ) :-- 



1 m ■ (e<°>) n -r<°) 

N L (i) . ( E (z)) N -R(z) 



where L^°\ R(°\ L^ T \ and R^ are implicitly func- 
tions of L s and R given by equation |2]). Assuming O 
corresponds to an extensive observable, we expect that 
limAr^oo £ (Ls, Rs) = (O) - that is, in the infinite limit 
the expected value per site converges to some number (O) 
which is independent of the boundaries. With this phys- 
ically reasonably assumption, we can reason about the 
structure of E^ and E^ to obtain an algorithm for 
computing (O). 

We first observe that since the maximal eigenvalue of 
E^ is 1 (due to the normalization of the site tensor, S), 
so must be the maximal eigenvalue of E^°\ since other- 
wise E N would be exponential in N . Furthermore, since 
£ N is linear in N for large N, the maximal eigenspace of 
must have a Jordan block structure-that is, there 
must be a matrix U = [u\ U2] with orthonormal columns 
Ui that provide a basis for this eigenspace such that 

1 



.4 



U.E<°)-W 



The matrix element a can be 



thought of as giving us the unnormalized expected value 
of O per site. In order to normalize it, we observe that 
as N — > co, 



?N 



1 L<°> • (uiu{ + u 2 u\ + Nauxul^j ■ R {0) 



N L {i) . (£<z>) 



N 



L (T) . (E(i)) N ■ RP) 

Since we expect £ N to be independent of the boundaries 
for large N, it must be the case that 



E&) -R^ 



N 



p\ L m.(^).rW , 



where 
u k ■ R 



-- u k ■ L 



o — 



Ei" u k. 



V") 



R% 



(i 1 ,i" ? 

that 



L%, and 



and 



are the projections of the eigenspace of E^ into a 
space applicable to E^ by dotting out the respective 
left and right boundaries of the matrix product opera- 
tor O; put another way, if we let V L := [v^ v^] and 
V R :— [v R v R ], then the limiting action of E^ in the 
(projected) eigenspace is given by B — 1 " L ■ r ' 1 ' ■ 1 ~ J!< 
[3 




V L. E (l).yi 

The matrix element (3 gives us the normalization 



constant, so that (O) := liniAr^oo £ = a//3. 

It remains to extract a and /3 from the respective ma- 
trices A and B. Our eigenvalue solver is not guaranteed 
to give us a basis of the eigenspace that makes A and 
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1. Compute the maximal two-dimensional eigenspace of 

— for example, by using ARPACK— and requesting a 
Schur basis rather than a (non-existant) eigenvector basis — 
and store the basis vectors in the columns of a matrix U. 

2. Express the action of in this space by computing the 
2 x 2 matrix A := U ■ E io) ■ UK 

3. Extract the unnormalized expected value from this matrix 
by computing a := \/tr (At • A) — 2. 

4. Project out the operator-dependent boundary conditions 
from the vectors in this eigenspace to obtain V L := L° -U = 

£i#'tfci',i",i'"),* and yR ~R° -u = £-R?'U(i'.i» *").*• 

5. Express the action of in this projected space by com- 
puting the 2 x 2 matrix B ~ V L ■ E {1) ■ V Rt . 

6. Extract the normalization constant by computing /3 := 
\/tr(Bt . B ). 

7. Obtain the (normalized) expected value by computing 
(O) := a/p. 
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FIG. 1: (color online) Finite state automata representation of 
a sum of exponentially decaying XX interactions. 



TABLE II: Algorithm to compute (O) for general matrix 
product operator O. 

B have the nice form above, but happily it is easy to 
see that a = y^r (A • At) - 2 and (3 = yjti (B -St), and 
these formulas are invariant under similarity transforms 
of A and B so they work independent of the basis we are 
given. The algorithm we have derived in the preceding 
paragraphs is summarized in Table [II] 



III. MATRIX PRODUCT FACTORIZATION OF 
EXPONENTIALLY DECAYING INTERACTIONS 

In the algorithms of Sec. [Til we assume that there 
is a matrix product representation of our Hamiltonian. 
It has previously been shown that there are matrix 
product factorizations of Hamiltonians with short-range 
interactions. 11 ' 12 In this section, we extend these results 
to show that there are also matrix product factorizations 
of Hamiltonians with long-range exponentially decaying 
interactions as well. 

So consider a Hamiltonian of an iV-site system experi- 
encing a sum of long-range exponentially decaying inter- 
actions, 

H:= J2 fe"n/£ J I® l ®X®I® r ®X®I^, 

i+l+r+l+j=N \ n ) 

where I is the identity operator and X is the Pauli opera- 
tor o~ x . We shall now employ the formalism in Ref. [l2| to 
obtain a factorization of this Hamiltonian. We start by 
changing our interpretation of the Hamiltonian: rather 
than thinking of it as a sum of tensor product terms 
with complex coefficients, we shall instead think of it as 
a function which maps arbitrary strings of X and I sym- 
bols to complex numbers. In our sum each term takes the 
form (J2 n u n (3 r n ) I® 1 <g> X <g> l® r <g> X ® 1® J with i,j,r>0, 



so our corresponding function maps strings of the form 
V X I r X V to Y^ n a nPn: and all other strings to zero. 

We shall now construct a finite state automaton which 
computes this function. A finite state automaton can 
be thought of as a machine which reads through an in- 
put string and changes its state in response to each in- 
put symbol according to a set of transition rules. The 
transitions are non-deterministic and weighted - that is, 
the automaton can take many transitions simultaneously, 
and at the end of the string it outputs a number corre- 
sponding to the sum of the product of the weights along 
each sequence of transitions that it took. The automa- 
ton starts on a designated initial state; if it does not end 
on a designated accept state or it encounters a symbol 
without an associated transition, then it outputs a zero 
(overriding other weights). 

The automaton which computes the function repre- 
senting our operator is illustrated in Fig. [TJ States are 
indicated by circles, and transition rules are indicated 
by arrows labeled with a symbol and a weight (one if 
not otherwise specified). An unconnected arrow desig- 
nates the initial state (1), and shading designates the 
accept state (2). A good way to think about what is 
going on is that terms are generated by each possible 
walk from state 1 to state 2. So for example, by tak- 
ing the path 1^1^3^3^3^2 we generate the term 
(ai/3 2 )I CS> X ® I <S> I <8> X; summing over all walks for 
strings of the form I X I 2 X we obtain the desired coeffi- 

cient J2n a nPl- 

From this automaton, we immediately obtain a matrix 
product operator representation. The initial and accept 
states give us the values for the left and right bound- 
aries: (L°), = Si : k and (R°) k = &i,k- The elements 
of the operator tensor, O^ , are given by the weight on 
the i — y j transition with the symbol corresponding to 
operator A (zero if no such transition exists); for exam- 
ple, we have that Of~ 3 = ai. To get a feeling for why 
this works, observe that a run of the automaton is equiv- 
alent to starting with a vector giving initial weights on 
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the states, multiplying this vector some number of times 
by a transition matrix, and then dotting the result with 
a vector that filters out all but the weights on the accept 
states; this procedure is exactly equivalent to the form of 
equation {T]). 

This process is easily extended to include terms with 
arbitrary spin coupling interactions, such as a x a Y , 
a Y a Y , g z ' a z , etc. Furthermore, one can combine a 
sum of several such interactions into a single automaton 
by having them all share the same starting and ending 
states. 



IV. RESULTS: HALDANE-SHASTRY MODEL 

Now we pull all of the ideas from the previous sections 
together and apply them to tackle the Haldane-Shatry 
model. 16,17 In the infinite limit this model is given by 
the Hamiltonian H. = &i ■ cfi+ r /r 2 , which features 

an anti-ferromagnetic dipole interaction which falls off 
with the square of the distance between sites. (Note that 
since this model features anti-ferromagnetic interactions, 
we need to work with blocks of two sites, as discussed 
at the end of section Hi Al ) Although H cannot be ex- 
pressed exactly as a matrix product operator, we can 
approximate it arbitrarily well by a sum of exponentially 
decaying interactions, H w £\ ct, • a i+r (£)„ a„/3^ _1 ) 
(with a n £ R and |/3 ra | < 1), which can be factored ex- 
actly using the technique in section ITTT1 Since there are 
three spin-coupling interactions, a x a x _ r , afcrY+ri an< ^ 
&fcrf +rl which we combine into a single automaton as 
discussed at the end of section HTT1 we obtain an automa- 
ton with a number of states equal to three times the num- 
ber of terms in the expansion, N, plus two more for the 
starting and ending states; this quantity gives us the size 
of the auxiliary dimension for the corresponding matrix 
product operator, c = 3N + 2. 

It remains to find the coefficients in this expansion. 
One approach is to numerically solve for the coefficients 
which minimize the sum of the squares of the difference 
between the approximation and the exact potential for 
distances up to some cut-off - that is, to find the mini- 
mizer of the function, 

N I'cutoff / 1 \ 2 

i=l r =l ^ ' 

where r cuto ff should be chosen to be just beyond the max- 
imum effective range of the approximation, since larger 
values of r cuto fj result in a longer running time for the 
minimization without resulting in a better fit. 

For our application of the algorithm, we used a nonlin- 
ear least-squares minimization routine from MINPACK to 
find coefficients for expansions with three, six, and nine 
terms; the resulting approximate potentials are plotted 
along side the exact potential in Fig. [3J The upper cut- 
off on r was set to 10000 because, as can be seen in Fig. 
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FIG. 2: (color online) Difference between the energy of the 
computed state and the energy of the exact ground state, 
plotted for each of the three exponential expansions that were 
employed. 

[21 this was just beyond the effective maximum range ob- 
tainable from a 9-term approximation. 

Given this approximate matrix product factorization, 
we applied the algorithm in Table [I] to compute a trans- 
lationally invariant matrix product state representation 
of the ground state for selected values of x^ employing 
each of the three-, six-, and nine-term expansions. The 
energy per site was computed using the algorithm in Ta- 
ble HH and compared to the exact value obtained from 
Ref. Qj]. The difference between these values (i.e., the 
residual) is plotted for each expansion as a function of 
X in Fig. [21 Note that the residuals for all three ex- 
pansions agree up to some point, and then diverge to 
different "floors" . This is because at first the small value 
of x is the dominating factor which limits the fidelity of 
the ground state, and then later as x becomes large the 
finite number of terms in the exponential approximation 
becomes the dominating factor. 

For the nine-term expansion, we also com- 
puted the two-point correlator - that is, C(r) := 
(a x ® J®( r - 1 ) (g) cr x \ - using the algorithm given in 
section Hi B I for computing the expected value of a local 
operator. The result for several values of x is p lotted 
in Fig. [3] against the exact value from Ref. 17. Note 
that our approximation gets good agreement up to 
some length, after which it becomes a constant. This 
is because the algorithm is attempting to approximate 
this correlator using a sum of decaying exponentials 
plus a constant term, analagous to how we used a 
sum of decaying exponentials to approximate the 1/r 2 
interactions; by increasing Xi we are increasing the 
number of terms available to track the correlator, which 
results in systematic improvement. 

V. CONCLUSIONS 

To summarize, in this paper we have presented an algo- 
rithm for computing the ground state of infinite ID sys- 
tems. This algorithm differs from the iTEBD algorithm 7 
in that it uses a variational approach instead of imag- 



7 




Distance between sites 

FIG. 3: (color online) (top) Decaying 3-, 6-, and 9-term ex- 
ponential approximations to 1/r 2 potential, (bottom) Two- 
point correlator for selected values of x using the the 9-term 
exponential approximation; a residual of the correlator for 
each value of \ m plotted simultaneously and labeled by e x . 

inary time evolution, and from the PWFRGi^ in that 
it considers an infinite system from the start. Further- 
more, since the algorithm itself employs matrix product 
operators, it has the important advantage of being capa- 
ble of modeling long range interactions, and in particular 



any interaction which can be approximated by a sum of 
decaying exponentials. In order to benchmark the algo- 
rithm, we have computed an infinite MPS for the ground 
state of the Haldane-Shastry model. The corresponding 
two-point correlators are in remarkable agreement with 
the exact solution up to distances above a thousand spins. 

In conclusion, our results indicate that this algorithm 
adds significantly to the existent tools to address ID 
many-body systems since it allows the properties of bulk- 
scale materials to be studied for realistic long-range po- 
tentials. Furthermore, it admits a natural extension to 
lattice systems in higher spatial dimensions, for which 
work is currently in progress. 

We note that upon completion of this work, we learned 
of simultaneous work on an equivalent algorithm by Ian 
McColloch in Ref. |l8|; his presentation includes a detailed 
comparison of the convergence of the iTEBlil and the 
variational approaches. 
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